function allCovs = GenerateCovariates_3column(saveBool, sequenceFiles, stimDur)
%
%
%   _______________________________________________
%   by Marcelo G Mattar (07/24/2012)
%   mattar@sas.upenn.edu


%% CHECK INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 3
    stimDur = 1.600;
end

if nargin < 2
    sequenceFiles = {...
        %{
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Experiment/Config/seqs/3targetsA.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Experiment/Config/seqs/3targetsB.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Experiment/Config/seqs/3targetsC.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Experiment/Config/seqs/4targetsA.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Experiment/Config/seqs/4targetsB.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Experiment/Config/seqs/5targetsA.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Experiment/Config/seqs/5targetsB.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Experiment/Config/seqs/5targetsC.txt' ...
        %}
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis/Covariate design/v1_0/Face1.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis/Covariate design/v1_0/Face2.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis/Covariate design/v1_0/Face3.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis/Covariate design/v1_0/Grating1.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis/Covariate design/v1_0/Grating2.txt' ...
        '/Users/marcelomattar/Dropbox/Marcelo/UPenn/Documents/Projects/AdaptID/Analysis/Covariate design/v1_0/Grating3.txt' ...
        };
else
    if ~iscell(sequenceFiles)
        sequenceFiles = cellstr(sequenceFiles);
    end
end

if nargin < 1
    saveBool = 0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% General
%expInfo.expType = 'fMRI'; % Experiment type. Enter "fMRI" or "Behavioal" to specify the directory where the data will be saved

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



numFiles = length(sequenceFiles);

% Initialize covariates
mainEffect = cell(length(sequenceFiles),1);
targetTrials = cell(length(sequenceFiles),1);
afterTargets = cell(length(sequenceFiles),1);
afterBlanks = cell(length(sequenceFiles),1);
carryOver1 = cell(length(sequenceFiles),1);


for seqnum = 1:numFiles
    
    % Load the sequence
    thisSeq = double(importdata(sequenceFiles{seqnum}));
    
    % mainEffect: non-blank trials have a value of 1; blank trials have a
    % value of 0
    mainEffect{seqnum} = ~(thisSeq==0);
    
    % Determine the trials following a blank block
    afterBlanks{seqnum} = zeros(length(thisSeq),1);
    afterBlanks{seqnum}(findEndOfBlankBlock(thisSeq)+1) = 1;
    afterBlanks{seqnum}(afterBlanks{seqnum} ~= 1) = -1;
    
    % Determine the target trials
    targetTrials{seqnum} = zeros(length(thisSeq),1);
    targetTrials{seqnum}(thisSeq==-1) = 1;
    targetTrials{seqnum}(targetTrials{seqnum} ~= 1) = -1;
    
    % Determine the trials following a target trial
    afterTargets{seqnum} = zeros(length(thisSeq),1);
    afterTargets{seqnum}(find(thisSeq==-1)+1) = 1;
    afterTargets{seqnum}(afterTargets{seqnum} ~= 1) = -1;
    
    % Calculate the difference between each stimulus and the previous one
    carryOver1{seqnum} = thisSeq - [thisSeq(end); thisSeq(1:(end-1))];
    carryOver1{seqnum}(carryOver1{seqnum} ~= 0) = 1;
    carryOver1{seqnum}(carryOver1{seqnum} ~= 1) = -1;
    
    
    
    % Remove any non-typical trial from the covariates
    removeTrials = (mainEffect{seqnum} == 0) + (afterBlanks{seqnum} == 1) + (targetTrials{seqnum} == 1) + (afterTargets{seqnum} == 1);
    
    afterBlanks{seqnum}(find(removeTrials)) = 0;
    afterBlanks{seqnum}(findEndOfBlankBlock(thisSeq)+1) = 1;
    
    targetTrials{seqnum}(find(removeTrials)) = 0;
    targetTrials{seqnum}(thisSeq==-1) = 1;
    
    afterTargets{seqnum}(find(removeTrials)) = 0;
    afterTargets{seqnum}(find(thisSeq==-1)+1) = 1;
    
    carryOver1{seqnum}(find(removeTrials)) = 0;
    
    
    
    % Mean center each covariate
    thisMean = sum(afterBlanks{seqnum}(afterBlanks{seqnum} ~= 0))/sum(afterBlanks{seqnum} ~= 0);
    afterBlanks{seqnum}(afterBlanks{seqnum} ~= 0) = afterBlanks{seqnum}(afterBlanks{seqnum} ~= 0) - thisMean;
    
    thisMean = sum(targetTrials{seqnum}(targetTrials{seqnum} ~= 0))/sum(targetTrials{seqnum} ~= 0);
    targetTrials{seqnum}(targetTrials{seqnum} ~= 0) = targetTrials{seqnum}(targetTrials{seqnum} ~= 0) - thisMean;
    
    thisMean = sum(afterTargets{seqnum}(afterTargets{seqnum} ~= 0))/sum(afterTargets{seqnum} ~= 0);
    afterTargets{seqnum}(afterTargets{seqnum} ~= 0) = afterTargets{seqnum}(afterTargets{seqnum} ~= 0) - thisMean;
    
    thisMean = sum(carryOver1{seqnum}(carryOver1{seqnum} ~= 0))/sum(carryOver1{seqnum} ~= 0);
    carryOver1{seqnum}(carryOver1{seqnum} ~= 0) = carryOver1{seqnum}(carryOver1{seqnum} ~= 0) - thisMean;
    
    
    
    if saveBool
        if ~exist('./covariates/', 'dir')
            mkdir('./covariates/');
        end
        i = max(strfind(sequenceFiles{seqnum}, '/'));
        j = max(strfind(sequenceFiles{seqnum}, '.'));
        sequenceName = sequenceFiles{seqnum}((i+1):(j-1));
        
        mainEffect_fileID = fopen(['./covariates/' sequenceName '_mainEffect.txt'],'w');
        targetTrials_fileID = fopen(['./covariates/' sequenceName '_targetTrials.txt'],'w');
        afterTargets_fileID = fopen(['./covariates/' sequenceName '_afterTargets.txt'],'w');
        afterBlanks_fileID = fopen(['./covariates/' sequenceName '_afterBlanks.txt'],'w');
        carryOver1_fileID = fopen(['./covariates/' sequenceName '_carryOver1.txt'],'w');
        
        for trial=1:(length(thisSeq)-1)
            fprintf(mainEffect_fileID,'%f\t%f\t%f\n',stimDur*(trial-1),stimDur,mainEffect{seqnum}(trial));
            fprintf(targetTrials_fileID,'%f\t%f\t%f\n',stimDur*(trial-1),stimDur,targetTrials{seqnum}(trial));
            fprintf(afterTargets_fileID,'%f\t%f\t%f\n',stimDur*(trial-1),stimDur,afterTargets{seqnum}(trial));
            fprintf(afterBlanks_fileID,'%f\t%f\t%f\n',stimDur*(trial-1),stimDur,afterBlanks{seqnum}(trial));
            fprintf(carryOver1_fileID,'%f\t%f\t%f\n',stimDur*(trial-1),stimDur,carryOver1{seqnum}(trial));
        end
        fprintf(mainEffect_fileID,'%f\t%f\t%f',stimDur*(length(thisSeq)-1),stimDur,mainEffect{seqnum}(length(thisSeq)));
        fprintf(targetTrials_fileID,'%f\t%f\t%f',stimDur*(length(thisSeq)-1),stimDur,targetTrials{seqnum}(length(thisSeq)));
        fprintf(afterTargets_fileID,'%f\t%f\t%f',stimDur*(length(thisSeq)-1),stimDur,afterTargets{seqnum}(length(thisSeq)));
        fprintf(afterBlanks_fileID,'%f\t%f\t%f',stimDur*(length(thisSeq)-1),stimDur,afterBlanks{seqnum}(length(thisSeq)));
        fprintf(carryOver1_fileID,'%f\t%f\t%f',stimDur*(length(thisSeq)-1),stimDur,carryOver1{seqnum}(length(thisSeq)));
        
        %{
        dlmwrite(['./covariates/' sequenceName '_mainEffect.txt'],mainEffect{seqnum});
        dlmwrite(['./covariates/' sequenceName '_targetTrials.txt'],targetTrials{seqnum});
        dlmwrite(['./covariates/' sequenceName '_afterTargets.txt'],afterTargets{seqnum});
        dlmwrite(['./covariates/' sequenceName '_afterBlanks.txt'],afterBlanks{seqnum});
        dlmwrite(['./covariates/' sequenceName '_carryOver1.txt'],carryOver1{seqnum});
        %}
    end
end

allCovs.mainEffect = mainEffect;
allCovs.targetTrials = targetTrials;
allCovs.afterTargets = afterTargets;
allCovs.afterBlanks = afterBlanks;
allCovs.carryOver1 = carryOver1;




function endOfBlankBlocks = findEndOfBlankBlock(sequence)

    startOfBlankBlocks = [];
    endOfBlankBlocks = [];
    
    if sequence(1) == 0; % If the first element is a blank
        startOfBlankBlocks = 1;
    end
    
    for i=1:(length(sequence)-1)
        if (sequence(i)~=0) && (sequence(i+1)==0)
            startOfBlankBlocks = [startOfBlankBlocks, (i+1)];
        end
        if (sequence(i)==0) && (sequence(i+1)~=0)
            endOfBlankBlocks = [endOfBlankBlocks, i];
        end
    end
